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Abstract 

In this paper we present a decentralized algorithm to estimate the eigenvalues of the Laplacian matrix that encodes the network 
topology of a multi-agent system. We consider network topologies modeled by undirected graphs. The basic idea is to provide 
a local interaction rule among agents so that their state trajectory is a linear combination of sinusoids oscillating only at frequen- 
cies function of the eigenvalues of the Laplacian matrix. In this way, the problem of decentralized estimation of the eigenvalues 
is mapped into a standard signal processing problem in which the unknowns are the finite number of frequencies at which the 
signal oscillates. 
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1 Introduction 

Nowadays the research community is investing more and 
more effort in designing coordination and estimation al- 
gorithms for networked multi-agent systems [1,4,8,16]. 
The network topology of a multi-agent system can be 
effectively described by means of a graph, where nodes 
represent agents and edges represent couplings between 
them [14]. The emergent behavior of such a system de- 
pends both on the interactions between agents and on 
the network topology. 

Algebraic graph theory [7] provides powerful tools to an- 
alyze a graph. As an example, the knowledge of the spec- 
trum of the Laplacian matrix associated to a graph can 
be used to estimate topological properties of an undi- 
rected graph, e.g., algebraic connectivity, average de- 
gree, diameter, spectral gap, connectivity measures [13]. 
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In the context of multi-agent systems, this knowledge 
may also provide powerful insights into the dynamics of 
a networked control system. As an example, as stated 
in [16], the algebraic connectivity, i.e., the second small- 
est eigenvalue, is a fundamental parameter to estimate 
the worst case convergence rate of consensus algorithms 
and, more in general, of multi-agent systems with local 
interactions described by the Laplacian matrix such as 
leader-follower networks [17]. 

Unfortunately, the spectrum of the Laplacian matrix is 
not readily computable in a distributed setting where 
the network topology is unknown. In order to overcome 
this limitation, we have designed a local interaction rule 
so that the resulting dynamical system oscillates only 
at frequencies corresponding to the eigenvalues of the 
Laplacian matrix that encodes the network topology. In 
this way, the problem of estimating the eigenvalues is 
mapped into a signal processing problem solvable inde- 
pendently by each agent in a decentralized fashion, us- 
ing tools from signal processing or system identification 
theory. 

Compared to the state of the art, discussed in detail in 
next section, our approach allows to estimate the full 
spectrum of the symmetric Laplacian matrix without 
the need to estimate all the corresponding eigenvectors. 
Moreover our approach provides an approximate esti- 
mation of the eigenvalues in finite time. 
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The contributions of this paper are the following: 

• We propose a novel local interaction rule to make the 
network state oscillate at frequencies corresponding to 
the eigenvalues of the Laplacian matrix, thus mapping 
the decentralized eigenvalue estimation problem into 
a standard signal processing problem. 

• We extend [5] by characterizing analytically the am- 
plitude and phase of the oscillations as function of the 
eigenvectors of the Laplacian and the initial condi- 
tions. 

• We propose an improvement with respect to [5] so that 
no component at null frequency exists in the evolu- 
tion of the agents' state. The removal of the DC com- 
ponent allows the straightforward application of fre- 
quency estimation algorithms such as the one in [3]. 

Related works 

In [21] Zavlanos et al. investigated the problem of how 
to coordinate a network of mobile robots with position- 
dependent topology so that the corresponding adjacency 
matrix has a given set of eigenvalues. This approach is 
based on artificial potentials, function of the inter-agent 
distances, that allow a gradient descent algorithm to 
make the network converge to a topology whose eigen- 
values are the desired ones. In this preliminary paper 
the authors consider the spectral moments related to the 
spectrum of the adjacency matrix to be known. 

In [6] Franceschelli et al. presented a necessary and suffi- 
cient condition to verify observability and controllability 
with of a leader- follower network with unknown topol- 
ogy of mobile vehicles is proposed based on the algorithm 
in [5] and its extension in this paper. 

In [20] Yang et al. proposed a technique for the esti- 
mation of the second smallest eigenvalue of a weighted 
Laplacian matrix based on the power iteration algo- 
rithm by the estimation of the corresponding eigenvec- 
tor. In addition, the authors discuss a decentralized con- 
trol algorithm to maximize the algebraic connectivity. 
The idea is to let agents move so that links are added or 
weights changed as two agents come closer. 

In [18] Sahai et al. presented an approach building on the 
idea of [5] for the application of clustering. The authors 
propose a local interaction rule formally equivalent to 
the wave equation discretized in time and space and show 
that the "wave equation method", can be used to cluster 
a graph by estimating the sign of the coefficients of the 
discrete Fourier transform corresponding to the second 
smallest eigenvalue. Furthermore, their approach is su- 
perior with respect to the convergence speed to the state 
of the art of eigenvector based clustering algorithms. 

Finally, in [9] Kempe et al. are interested in computing in 
a decentralized way an approximation of the first k eigen- 
vectors of a symmetric matrix that encodes the network 



topology. Their algorithm takes inspiration from the or- 
thogonal iteration algorithm and assumes that the net- 
work topology is unknown to the nodes. This algorithm 
can be adapted to our objective, i.e., the distributed es- 
timation of the eigenvalues of the Laplacian matrix, by 
introducing a distributed technique for the computation 
of the Rayleigh quotient. 

2 Online Spectrum Estimation 

Let us consider the interactions of a network of agents 
described by an undirected graph Q = {V,£}, where 
V = {1, . . . , n} is the set of agents and £ C V x V is the 
set of edges: an edge eij exists between agents i and j if 
agent i interacts with agent j. 

Let Mi define the neighborhood of agent i, namely the 
set of indices of the agents connected by an edge with 
agent i. In particular, \Mi\ = Aj where Aj is called de- 
gree of agent i. Let C be the Laplacian matrix of graph 
Q, it is a n x n matrix the elements of which are Zy = Aj 
if z = j, lij = — 1 if j G Mi and otherwise. The Lapla- 
cian matrix C of an undirected graph is symmetric by 
construction and thus all its eigenvalues are real. Fur- 
thermore, for a connected graph, it has one null struc- 
tural eigenvalue with corresponding eigenvector equal to 
the vector of ones 1„ of appropriate dimensions, thus 
C 1„ = 0„. In addition, according to the Gershgorin disc 
theorem, a symmetric Laplacian has all its eigenvalues 
located within [0, 2 A max ], where A max is the maximum 
degree between the agents in the graph. 

We now present a decentralized algorithm to estimate 
the eigenvalues of the Laplacian matrix. The algorithm 
requires each agent i to store two variables Xi, Z{ £ M 
and apply a local update rule upon receiving the values 
of the equivalent variables from its neighbors. 

Algorithm 1 (Online Spectrum Estimation) 

(1) Each agent sets t = and chooses an initial condi- 
tion uniformly at random Xi(0), Zj(0) € {—1, 1}- 

(2) Each agent simulates the following local interaction 
rule with its neighbors Mi(t) 

Xi (t) = Zi (t) + ( Zi (t) - Zj (t)) , 

ii(t) = -Xi(t) - EjeJVi (^W - x i( t )) ■ 

(3) In a time window of length T , agent i estimates the 
frequencies of the sinusoids of which signal Xi(t) is 
composed. 

(4 ) The values of the frequencies estimated correspond 
to the eigenvalues of the Laplacian matrix C shifted 
by 1 and are given as output. M 

Note that Step 3 can be solved by several methods of 
signal processing or system identification. In particular, 
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the required value of T depends on the chosen algorithm. 
In this paper, as discussed in Section 3, we exploit the 
method presented in [3] to implement Algorithm 1. 

The behavior of the network when all the agents update 
their state according to eq. (1) can be described as fol- 
lows 
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where f is the frequency domain variable. In addition, 
the coefficients ajj and bj t i are given by 



- For Ai = 



ai,i = vi(i)vfx(0) 
bi.i = vi(i)viz(0) 



(3) 



where / is the n x n identity matrix and 0„ x „ is the 
null n x n matrix. Note that for any network topology 
A is skew symmetric, i.e., A T = —A. In the following 
theorem, we prove that the eigenvalues of A can be an- 
alytically derived from the eigenvalues of the Laplacian 
matrix C and they are all structurally purely imaginary. 

Lemma 1 Let Q be an undirected graph with Laplacian 
C. Let matrix A be defined as in eq. (2) . To any eigenvalue 
Xc of C it corresponds a couple of complex and conjugate 
eigenvalues X A , X A of A, that is: 

\A= jQ- + *£.), Aa = - ? (1 + A £ ), 

while the corresponding eigenvectors v\ A are function of 
the eigenvectors v\ c of C 



For Xj > 
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Proof: See Appendix B. 



(4) 



□ 



The above theorem states the key result of this paper. In 
fact, it implies that each agent can independently solve 
the problem estimating the eigenvalues by estimating 
the frequencies at which its own state variable Xi (t) os- 
cillates. 
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Proof: See Appendix A. 
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By Lemma 1 it follows that the state of each agent has an 
oscillatory trajectory. Furthermore, as detailed by Theo- 
rem 2, this trajectory is a linear combination of sinusoids 
oscillating only at frequencies function of the eigenval- 
ues of the matrix Laplacian. In the following we assume 
the Laplacian to have m distinct eigenvalues labeled as 
follows: = Ai < A 2 < • • • < X m . 

Theorem 2 Let us consider a system described by eq. (2) 
relative to a network whose graph Q is connected. Let 
x(0) = xq and z(0) = zq be the state initial conditions. 
LetS(-) be the Dirac's delta function. Let Xj be an eigen- 
value of the Laplacian matrix C of graph Q with alge- 
braic multiplicity i>j and let m be the number of distinct 
eigenvalues. Let V\ be the unitary norm eigenvector cor- 
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responding to Ai = 0, and v 
unitary norm eigenvectors associated to Xj > 0. The 
module of the Fourier transform of the i-th state compo- 
nents Xi(t) and Zi(t), i = 1, . . . , n, can be written as 



\Fix i (t)]\ = \Xi(f)\ = J2^r 6 (f ± 
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Remark 1 Few important remarks are now in order: 

- The value ofXi(t) can be seen as the output of the i-th 
agent. If the system is not observable from the output 
Xi(t) then some coefficients <Xj are null and thus the 
corresponding mode cannot be detected by agent i. 

- For each agent, the amplitude of the sinusoid oscillat- 
ing at uj — Ai = 1 corresponds to the instantaneous 
average of the state variables. M 

The observability and controllability of a system, the 
dynamics of which are described by the Laplacian ma- 
trix, have been studied in [11,12] from a graph theoret- 
ical point of view. In the following theorem we link the 
ability to estimate all the eigenvalues of the Laplacian 
to its observability property. 

Theorem 3 Let A be the matrix describing the group 
dynamics as in (2). Let C be a k x n output matrix, with 
k G N. Let 
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Let O a = 0(A, C) and O c = 0(C, C) be the observabil- 
ity matrices built with the corresponding matrices. Then: 

Rank (O a ) = 2Rank (O c ) ■ 
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Proof: See Appendix C. 



□ 



We now state the main result of this paper that proves 
the correctness of Algorithm 1 . 

Theorem 4 Consider a connected network Q ofn agents 
that executes Algorithm 1. Let the initial conditions of 
system (2) be not orthogonal to any eigenvector of matrix 
C. Let C = [0 . . . , 1, . . . , 0] be zero everywhere except for 
the i-th unitary element, with i £ V . If O (C, C) is full 
rank, then agent i can estimate all the eigenvalues of the 
Laplacian matrix. 

Proof: Due to Lemma 1 all the eigenvalues of system (2) 
are purely imaginary and correspond to the eigenvalues 
of the Laplacian matrix shifted by one. Furthermore, if 
the initial conditions are not orthogonal to all the eigen- 
vectors of C and O (C, C) is full rank as discussed in The- 
orem 3, then all the sinusoids corresponding to the sys- 
tem modes have coefficients strictly greater than zero. 
Thus by applying a frequency estimation algorithm to 
the signal for instance the one in [3], agent i can 

estimate the full spectrum of the Laplacian matrix by 
only observing its own state evolution. □ 

It is relevant to point out that even if the system is not 
observable from a single agent perspective, it will always 
be observable if matrix C is the identity matrix, i.e., if we 
consider all the information that agents locally retrieve. 





Fig. 1. Topology variation with respect to time for a network 
composed of 5 agents. 
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Fig. 2. Spectrogram of the time varying topology shown in 
Fig. 1 computed by the ith agent with respect to the output 
associated to its state variable Xi(t). 



3 Numerical implementation of the approach 

The system in eq. (2) is a marginally stable linear sys- 
tem since all its eigenvalues lie on the imaginary axis. 
The stability of a system with eigenvalues exactly on the 
imaginary axis is not considered to be robust because 
even the slightest parameter uncertainties may render 
the system unstable. In our case there is no parameter 
uncertainty because system (2) is based on the Lapla- 
cian matrix the elements of which depend only on the 
existence of links between the agents. Thus, for any net- 
work topology system (2) can not be stable or unstable 
but only marginally stable. Furthermore we point out 
that since no sensing/measurement is involved, no noise 
is generated from the application of the local interaction 
rule. 

In this paper, we implemented the approach in [3] to 
estimate the frequencies at which the signal oscillates. 
Furthermore we performed a spectral analysis by means 
of the Discrete Fourier Transform (DFT) . 

3.1 Approximate Frequencies Estimation Method 

The problem of estimating the frequencies of a signal 
such as 
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Fig. 3. Spectrum of the time varying topology shown in Fig. 1 
computed by the ith agent at time t = 4, t = 12, t = 20. 
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y(t) = Aisin(uit + (f>i) 



(5) 



Table 1 

Comparison between the actual spectrum of the Laplacian, 
denoted by Xc, and the estimated spectrum, denoted by Ac 
of the time varying topology shown in Fig. 1. 
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Fig. 4. Eigenvalues estimated every T = 1 units of time 
by the approximate frequency estimation method in Sub- 
section 3.1. Note that the frequencies are shifted by 1 with 
respect to the eigenvalues of C in Table 1. 

has been extensively studied in control theory via off- 
line methods based on Fourier analysis tools and on- 
line methods [10,15]. In this paper, the approximate 
frequency estimation algorithm in [3] has been imple- 
mented. It allows to estimate the frequencies of a sig- 
nal of the form in eq.(5) by assuming an upper bound 
nMax to the number of existing frequencies to be avail- 
able. The input of the algorithm is the sampling time 
T s , an upper bound nMax of the number of expected fre- 
quencies, and a measure of the approximation error S e . 
Furthermore, the length of the time window considered 
must be greater than the largest period of the sinusoid 
with the smallest frequency. 

The output of the algorithm is the number n of esti- 
mated frequencies and their values {1 + Xc} and a flag f 
the value of which is true if the percentage error when 
reconstructing the signal with the estimated coefficients 
fits the threshold S e , false otherwise. Note that, a great 
advantage of this algorithm is that the estimation can be 
worked out in finite time. If an observer with asymptotic 
convergence is required, the output of the algorithm de- 
scribed above could be used as input for the algorithm 
proposed in [2]. 

3.2 Simulations with switching topology 

In order to corroborate the mathematical results, simu- 
lations have been carried out by exploiting the 4th Order 
Runge-Kutta Method (RK4) to simulate the system (2). 
Regarding the signal processing, let us recall that this 
can be always carried out locally by each agent in spite 
of the particular technique adopted. 

In the simulation, a network of agents the topology of 
which changes over time is considered. In detail, Fig. 1- 
a) depicts the network topology at the time interval t £ 
[0,6.4), Fig. 1-b) describes the network topology at the 
time interval t £ [6.4, 12.9) and Fig. 1-c) describes the 
network topology at the time interval t <E [12.9, 20]. Each 
agent is running the interaction rule described in eq. (1). 

Fig. 2 shows the spectrogram of the time varying topol- 
ogy computed by the agent i with respect to its state 



variable Xi(t). The spectrogram was computed by this 
agent with f s = . The x and y axes of the spectro- 

gram represent respectively the time step and angular 
frequency, while the color of the spectral line describes 
the amplitude of the frequency peaks, i.e., white means 
zero amplitude while black means an amplitude greater 
than 0.1. 

Fig. 3 shows a section of the spectrogram at different 
time steps, namely t = {4, 12, 20}, representing the 
spectrum of the three network topologies taken into ac- 
count. 

To this example, we also applied the method in Sec- 
tion 3.1 to estimate the frequency of the sinusoids and 
thus the eigenvalues of the Laplacian. The comparison 
between the eigenvalues of the Laplacian matrix of the 
time-varying network topology in Fig. 1 and the esti- 
mated eigenvalues in Fig. 4 using the approximate fre- 
quency estimation method in Subsection 3.1 is shown in 
Table 1. This method was implemented choosing as sam- 
100 



pling frequency f s 



2tt 



that is twenty times the max- 



imum expected frequency in the signal. This frequency 
corresponds to the largest eigenvalue of the Laplacian 
matrix plus one. The length of the time window used 
to computed each estimation is T = 1 sec, which is the 
period of the slowest sinusoid. 

3.3 Computational Cost 

To study the computational cost of Algorithm 1 we adopt 
the metrics proposed in [9], i.e., we count the number 
of communication rounds required to obtain an estima- 
tion with a certain accuracy. In this view, two important 
aspects must be considered: (i) the proposed algorithm 
consists in a local interaction rule that is supposed to be 
applied continuously; (ii) the required signal processing 
is carried out locally by an agent and several techniques 
can be adopted. Therefore, the computational cost anal- 
ysis consists in the study of the cost of simulating the lo- 
cal interaction rule and the cost of the signal processing. 
Since the study of the computational cost for the sig- 
nal processing required to estimate a discrete number of 
frequencies contained in a signal is not the scope of this 
paper we focus our attention to the number of communi- 
cation rounds required for the discrete time simulation 
of system (2) by the agents to collect enough data for the 
consecutive signal processing. In particular, the simula- 
tion time needed to collect a sufficient amount of data 
must be greater than the largest period of the sinusoid 
with the smallest frequency, T m ; n . By considering an ac- 
curate numerical simulation method such as the fourth 
order Runge-Kutta method, 4 messages have to be ex- 
changed between each agent and any of its neighbors to 
compute a sample of the state trajectory. It follows that 
for each agent the rounds of communication required to 
collect a sufficient amount of data can be bounded by 
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*± Ll max ■T min -f s , with A max the maximum degree in the 
network and f s the chosen sampling frequency that has 
do be at least greater than twice the largest eigenvalue of 
the Laplacian, thus greater than 2A max , to avoid alias- 
ing issues. 

4 Conclusions 

In this paper a decentralized algorithm to estimate the 
Laplacian spectrum of an undirected graph has been pro- 
posed. Each agent interacts with its neighbors so that 
its state oscillates at the frequencies corresponding to 
the eigenvalues of the Laplacian matrix that encodes the 
network topology. Therefore, the problem of estimating 
the eigenvalues has been reduced to a simple and widely 
studied problem of signal processing which involves the 
estimation of the discrete number of frequencies at which 
the generated signal is oscillating. A theoretical analysis 
of the proposed technique along with numerical simula- 
tions has been provided. 
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A Proof of Lemma 1 

By definition, the eigenvalues of A are the solutions of 
det(A - XI) = det 



-XI I + C 
-I-C -XI 



= 0. 



Since A is a block matrix whose blocks commute [19], 
then 



det{A - XI) = det ( X 2 I -{I + C) 



(A.l) 



Now, denote Xc the generic eigenvalue of the matrix 
Laplacian, it is det(C — Xcl) — 0. By adding and sub- 
tracting the identity matrix and exploiting the fact that 
the eigenvalues of the square of a matrix are squared, we 
obtain 

det((I + C) 2 - {l + X c ) 2 I)=0. 



Thus, by (A.l), it is 

x 2 = -(i + x c y 



X = ±j(l + X c ). 
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Therefore, to each eigenvalue of the Laplacian matrix C 
it corresponds two imaginary and conjugate eigenvalues 
of matrix A denoted by A_4 and Xa , so that 

\4=j(l + A £ ), and Xa = —j (1 + Xc) , 

thus proving the first statement. Now, by definition, the 
eigenvectors of A relative to Xa are solutions of 



I. Thus 
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for which a possible solution is va 
The same argument holds for the conjugate eigenvalue 
Xa = —j (1 + Xc) f° r which a possible solution is 
"A = i v L -3' V \J ■ 

B Proof of Theorem 2 

When referring to the eigenvalues and eigenvectors of C, 
Xc.j and v\ c , for j = 1, . . . , n, we drop the subscripts C 
and Xc , respectively, and refer to them as Xj and Vj for 
j = 1, . . . , n. By Lemma 1 to each Laplacian eigenvalue 
Xj it corresponds a couple of pure imaginary eigenvalues 
of A equal to 

X A ,X A =±j(l + Xj). 

Therefore, the state trajectory of each agent is a 
linear combination of sinusoids whose amplitudes and 
phase shifts are function of the initial conditions and of 
the graph topology. 

Now, we compute the coefficients of the module of the 
Fourier transform of Xi(t). Since A is skew symmetric, 
and any skew symmetric matrix is a normal matrix, 
i.e. A* A = AA* , thanks to the Spectral Theorem it is 
always diagonalizable through a unitary matrisQ] and 
all the eigenvalues have geometric multiplicity equal to 
their algebraic multiplicity (or equivalently, unitary in- 
dex). Thus A can be decomposed as A = VDV* , where 
D is a diagonal matrix whose elements are arranged as 
L> = diag{jAi, jX n , -3X1, -jX n }, and V 

is a complex matrix whose columns are the eigenvectors 
of A. Furthermore, applying Lemma 1, matrix V is re- 
arranged to match the disposition of the eigenvalues of 
D as follows 



V = 



vi v 2 . . 
jvi jv 2 . . 



jv n -JVx 



I '2 

-JV2 



-JVn 



In the following it is assumed that Vj, j = 1, . . . ,n, are 
normalized eigenvectors such that ||uj|| = 1 and VV* = 



1 A unitary matrix U is a complex matrix such that U*U = 
UU* — I, where U* is the complex conjugate of U. 



■ T 



1, a G 



By simple manipulations we find a = The state tra- 
jectories of the system are captured by the matrix expo- 
nential which in our case takes the form e At = Ve Dt V* . 
It follows that the state trajectory of Xi(t) and Zi(t) have 
the following form 

n 

Xi(t) = ]T [vj(i) (cos((l + Xj) t)vjx(0) 
3=1 

+ sva((l + Xj)t)vjz(0))] , 

n 

*(*) = E M l ) (-sin((l + Xj) t)vjx(0) 
j'=i 

+ cos((l + Xj)t)vjz{0))] . 

Thus, according to the notation of Theorem 2, by sim- 
ple manipulations we can obtain the expression for the 
coefficients a^j and bjj associated to the eigenvalue Aj 
for the i-th agent stated in eq. (3) and eq. (4), proving 
the theorem. 

C Proof of Theorem 3 

Let O a = 0(A, C), O i+c = 0(1 + C,C) and O c = 
0(£, C) be the observability matrices of the correspond- 
ing matrices. It can be shown by row permutation that: 



Rank ( O a ) = 



Rank 



V 



Oi+c 

o 1+c (i + cy 

Onxn 
Onxn 



Onxn 
Onxn 

o I+c 

I+c (I + C) n 



hence it also holds that 



C/+£ Onxn 
Onxn C/+£ 



Rank (0.4 ) = Rank 
from which it follows that Rank (Oa) = 2Rank (Oi+c)- 



Finally, by noticing that the eigenvalues of the matrices 
Oi+c and Oc are related as Xi+c = 1 + Xc and share 
the same eigenvectors, from the PBH observability test 
it follows that 



thus 



Rank (Oc) = Rank (O i+c ) , 
Rank (O a ) = 2Rank (O c ) ■ 
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